Low rank perturbations and the spectral statistics of pseudointegrable billiards 



Thomas Gorin 

Theoretische Quantendynamik, Universitat Freiburg, Hermann- Herder-Str 3, D-79104 Freiburg 



Jan Wiersig 

Institut fur Theoretische Physik, Universitat Bremen, Kufsteiner Str.. 

(Dated: February 8, 2008) 



D-28359 Bremen, Germany 



o 
o 

(N 



00 

Q 

U 
d 



(N 
> 

m 
o 
m 
o 
m 
O 



> 
X 



We present an efficient method to solve Schrodinger's equation for perturbations of low rank. In 
particular, the method allows to calculate the level counting function with very little numerical effort. 
To illustrate the power of the method, we calculate the number variance for two pseudointegrable 
quantum billiards: the barrier billiard and the right triangle billiard (smallest angle 7r/5). In this 
way, we obtain precise estimates for the level compressibility in the semiclassical (high energy) limit. 
In both cases, our results confirm recent theoretical predictions, based on periodic orbit summation. 
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Consider the stationary Schrodinger equation for a 
quantum system, which can be described by a Hamil- 
tonian H = Ho + W with the following properties: The 
operators H and Ho have discrete spectra, the eigenbasis 
of Hq is known, and the perturbation W is non-negative 
(or non-positive) and of low rank. Then, our method al- 
lows to evaluate the level counting function at arbitrary 
energies, by solving an eigenvalue problem of the dimen- 
sion which is equal to the rank of the perturbation. At 
first sight, the requirements seem rather restrictive, but 
in fact, there are a number of problems which can com- 
ply with them. For example, one may think of a few 
particle system with some kind of short range residual 
interaction. 

To illustrate how the method works, we choose two 
examples of a different type, the so called "barrier bil- 
liard" , and the right triangle billiard with smallest 
angle 7r/5. Both examples are twodimensional, pseudoin- 
tegrable polygon billiards 0, Q . While pseudointegrable 
systems have enough constants of motion to assure local 
integrability, singularities in the Hamiltonian flow allow 
invariant surfaces with genus larger than one. 

In the spirit of quantum-classical correspondence, 
there have been numerous efforts to study the impli- 
cations of pseudointegrability on the quantum spec- 
trum 0, H H H Q, @. In Refs. d 13 it is conjec- 
tured that the statistical properties of pseudointegrable 
billiards is intermediate between Poisson statistics, typ- 
ically related to integrable systems and the eigen- 
value statistics of the Gaussian orthogonal ensemble [llj , 
related to fully chaotic (time reversal invariant) sys- 
tems. The so called "intermediate statistics" has also 
been found in disordered, mesoscopic systems at the 
metal-insulator transition [l^l j for systems with interact- 
ing electrons |13| , and for incommensurate double- walled 
carbon nanotubes If 4j . 

A suitable measure for intermediate statistics is the 



level compressibility x = liniL^co £ 2 (L)/L, where £ 2 (£) 
is the number variance for energy intervals of length 
L (measured in units of the average level spacing) . Note 
that x coincides with the value of the spectral two-point 
formfactor in the limit of small times. Recently, analyt- 
ical results for x became available for a certain class of 
right triangle billiards (including our example) [T^|. as 
well as for the barrier billiard |8||. However, the conver- 
gence to the semiclassical limit is so slow, that numerical 
studies could not confirm those results with the desirable 
clarity. For example in the case of the barrier billiard 
one obtains x ~ 0.34 in the region between level number 
400 000 and 420 000, while the semiclassical prediction is 
X = V2 0- 

As we will show below, our method is almost ideally 
suited for the calculation of the number variance for 
large values for L. Therefore, we are able to calculate 
the level compressibility at much higher energies and 
with better statistics. In the case of the barrier billiard, 
for example, we calculate number variances in the region 
of absolute level number Af > 1.6 * 10 7 within an energy 
range which contains about 10 5 levels. 

Here we can only sketch the general method. A de- 
tailed presentation will be published elsewhere, and an 
early implementation can be found in Ref. pj. As- 
sume, that the Hamiltonian H can be approximated by 
a projection on an appropriately chosen iV-dimensional 
Hilbertspace, N < oo. Furthermore, assume that its ma- 
trix representation is of the following form: 



H = Ho 



yyt 



(i) 
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Here V is a iVxM-matrix with mutually orthogonal col- 
umn vectors, M < N, and r\ is a positive parameter. 
Such a representation can be obtained, for instance, by 
diagonalizing a non-negative perturbation of small rank 
M . Then, in general, the spectrum of H consists of a 
trivial component So, contained in the spectrum of Ho, 
and a non-trivial component Si, which is disjoint. The 
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eigenvalues in Si are roots of the secular equation 



det[l - r)K(E)} = , K{E) = V f 



1 



E-Hr 



V, (2) 



where the matrix K(E) is of dimension M 0. Eq. can 
be considered as an eigenvalue equation K(E)x = 5 x, 
where one of the eigenvalues must be equal to rj . 

Differentiating K{E) with respect to the energy gives 
a negative definite matrix. This shows that the eigen- 
values of K{E) are monotonously decreasing with en- 
ergy. The eigenvalues of Hq coincide with the positions 
of the poles of K(E). Therefore, each time the energy 
moves across a pole, a new eigenvalue of K(E) appears 
at +oo. The eigenvalue decreases with energy, until it 
reaches the value 1/77. At this point the secular equa- 
tion J2| has a root. Beyond this point, the eigenvalue 
continues to decrease, until it disappears at —00. This 
behavior gives rise to the following sum rule: Let no(ni) 
denote the number of eigenvalues of K(E), larger than 
I/77, at E = Eo(Ei), and let N p denote the number of 
poles and N r the number of roots of det[K(E)] in the 
interval (Eq,Ei). Then it holds: 



n + N p -N r 



ni 



(3) 



At the one hand, this relation can be used to bracket the 
eigenvalues in Si , before using a standard root searching 
algorithm to obtain the desired accuracy. In this case, 
Eq. © assures, that no eigenvalues are overlooked. 
On the other hand, in order to study the long range 
correlations, it can be sufficient to calculate the level 
counting function at the points of an equi-distant grid 
(in units of the average level spacing) grid in the energy 
region of interest. If L st denotes the stepsize, one can 
obtain the number variance at integer multiples of L st 
from this data. 

To calculate the spectrum of the barrier billiard 8], 
we choose as Ho the Hamiltonian of a rectangle billiard 
with sides of length a and 6. With the origin of a Carte- 
sian coordinate system fixed at one of the corners, the 
normalized eigenfunctions are 

f mn (i,i/) = 2(ab)~ 1 ' 2 sm(irmx/a) sin(7rm/ /b) , (4) 

while e mn = [(m/a) 2 + (n/6) 2 ]7r 2 /2 are the correspond- 
ing eigenvalues. We use units, in which the mass and 
Planck's constant h are both equal to one. The pertur- 
bation consists of an additional boundary segment inside 
the billiard, connecting the points (ao,0) and (ao,c). It 
is modeled by a potential well with delta shaped profile: 

H = H + r)W, W = (a/2) 5(x - ao) 0(c - y), (5) 

where 5(x) is the usual 5-function, and 9(y) is the unit 
step function. The prefactor a/2 is included for conve- 
nience. As rj increases from to 00, the spectrum of H 
changes from the spectrum of the rectangle billiard to 
the spectrum of the barrier billiard. 



In what follows we set ao = a/2, c = 6/2, where 
a = 27r 3 / 2 /3 and b = 6/71" 1 / 2 . In this way, we obtain 
the same spectrum as in Ref. In this case the triv- 
ial component So of the spectrum corresponds to states 
^mn with a node line at x — a/2. In order to obtain the 
non-trivial component alone, we require the eigenstates 
to be reflection symmetric with respect to that line. In 
other words, we consider the spectrum of a new rect- 
angle billiard with sides a/2 and b, which has Dirichlet 
boundaries everywhere, except for the boundary segment 
between the points (a/2,6/2) and (a/2,6) where it has a 
von Neumann boundary. 

To obtain the decomposition W — VV T , we calculate 
all eigenvectors of W which correspond to nonzero eigen- 
values. This is rather simple because W is separable in 
the x- and y-modes of the eigenfunctions of Ho : 

V™ = s(m) s(m) = sin(7rm/2) (6) 

V2a (_)«+« 
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(7) 



For the matrix elements of K{E) we get: 
v (a) V (J3) 

K ap (E) = Y, Y ^T^=Y, A( n ) W 



run mn 



s{m) 2 

E Em. 7 



(8) 

In the limit m — > 00, we can evaluate the sum over m 
analytically: 

K aff (E) = ^-J2 A n )A i P)G n(E), (9) 
n 

with G n (E) = — taxi(-Kz n /2) / z n . Here, z n is the effec- 
tive quantum number for Ho at given energy E, i.e. 
(z n ir/a) 2 + (nn/b) 2 — 2E. Note that G n (E) remains 
real, even for imaginary z n . Introducing the orthogonal 
matrix A nm = \/2A^ l _ 1 we may write: 



K{E) 



An 



(A G oM A 1 + G cvcn ) , 



(10) 



with G odd = diag[G 2 n-i(£0], and G cvcn = diag[G 2 „(£)]. 
Multiplying Eq. I|1U|) from left and/or right by A T and/or 
A, one can construct additional variants: L = A T K, 
K' = A T K A, and M = K A. Any one can be used in 
the secular equation to find the eigenvalues of the 
barrier billiard. However, the matrices K(E) and K'(E) 
are real and symmetric, and hence easy to diagonalizc. 
For the analysis below, we use K'(E), because it turned 
out that the numerical procedure converges faster in this 
case. In Ref. |8j the non-symmetric matrix M(E) has 
been used, instead. 

In what follows, we aim at a precise, numerical esti- 
mate for the level compressibility, which can be compared 
to the analytical results obtained in To this end, 

we compute the level counting function on a finite grid 
No, itotj -^st of consecutive intervals of length L st , start- 
ing at No and ending at No + L to t ■ The grid is defined on 
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(A) (B) (C) (D) 
M) = 8 * 10 5 2 * 10 6 8 * 10 6 1.6 * 10 7 
Ltot = 10 5 1.4 * 10 5 10 5 10 5 



TABLE I: Data sets for the barrier billiard. 

the unfolded energy axis, where the average level spacing 
is equal to one. It is mapped onto the physical energy 
axis, by inverting the Weyl law |l6j |: 

We count the number of levels in each interval, using the 
sum rule, Eq. |JSJ. The data sets, listed in Tab. [I] are 
produced in this way. Finally, they are used, to calculate 
the number variance E 2 (L) at integer multiples of L st = 
10. To get an idea of the efficiency of the method, note 
that it took about twenty days on a PC running at 1.7 
GHz, to produce the data set (D). The dimension of the 
matrices, to be diagonalizcd, was approximately 3700. 

According to Ref. |ljj, the number variance S 2 (L) sat- 
urates at L « L maXl which is related to the inverse period 
of the shortest (non-diffractive) periodic orbit. For our 
barrier billiard 

L max = y/2n Wo b/a = y/18 N /n . (12) 

At large L>1, the ratio Y?(V)jL is constant at first, but 
it starts to decrease, when L approaches L max . There- 
fore, we may expect that a plot: Y?{L)/L versus L/L max 
depends only weakly on L max , i.e. the energy range jVo- 
We call this function the scaled number variance: 

g(x) = ¥. 2 (L max x)/(L max x) . (13) 

In theory it is expected that in a given energy region, the 
ratio Y?{L)jL approaches its limit value quite quickly; 
for instance, L sw 10 seems to be sufficient, to obtain 
accurate results. The semiclassical limit Afo — * oo is the 
really difficult one. In this limit the level compressibility 
must be estimated by extrapolation: x — bur^o d( x )- 

Fig. n shows the scaled number variance for the data 
sets given in Tab.Q] For the error bars, we estima ted that 
the relative error is approximately equal to 1.6 sjL/ Ltot • 
In order to check this, we computed the variance of the 
distribution of different partial averages, fitting them 
with a normal distribution. Note that L to t/L gives the 
number of independent level counts, in the energy range 
considered. Fig. ^ explains the discrepancy between the 
numerical estimate for the level compressibility obtained 
in Ref. Q , and the theoretical expectation. In the region 
x > 0.01 the scaled number variance looks perfectly lin- 
ear. Hence, in the absence of data points with x < 0.01, 
one is lead to assume that the linear behavior continues 
down to the point x — 0. However, for x < 0.01, the slope 
changes drastically, and a second linear regime appears. 
Then indeed, the scaled number variance g(x) approaches 
the theoretical prediction x = 1/2 as x -> 0. In order 
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FIG. 1: The scaled number variance E 2 (L)/L versus x = 
L/Lmax, for the four data sets given in Tab. [I] For clarity, 
only some of the data points are plotted. The thick solid 
line gives the fit function g(x) with parameters as given in 
Eq. flg) . 

to put our findings on a quantitative basis, consider the 
following phenomenological parametrization: 

g(x) =a,2 — a\X + (ao - a 2 ) exp(-a 3 x) . (14) 

Its form is such that ao = <?(0) gives the best estimate 
for the level compressibility in the semiclassical limit. At 
the same time, both linear regimes are reproduced. By 
consequence, a-i would be the estimate for the level com- 
pressibility, in the absence of data points with x < 0.01. 
For the fit, all data sets in Tab. [i] are taken into account, 
excluding only those data points, for which x > 0.2. 
Beyond this point (which is outside the interval shown 
in Fig. ^) the parametrization l|14f> breaks down. Us- 
ing the nonlinear least squares Marquardt-Levenberg al- 
gorithm 0| , we obtain the following estimates for the 
fit-parameters: 

a = 0.5008(91) a x = 1.077(18) 

a 2 = 0.3866(17) a 3 = 364(39) (15) 

with a reduced chi-square value of xit// ~ 0-^4 (the 
number of degrees of freedom is / = 305). 

With ao, we have finally obtained a precise numerical 
estimate for the level compressibility in the semiclassical 
limit. It agrees within a relative error of roughly two 
percents with the theoretical result x = 1/2- In addition, 
the moderate value for x\t IS quite remarkable. It gives 
some support to the assumption that within the statis- 
tical error, the scaled number variance is independent of 
the energy region (as long as x < 0.2). 

In what follows, we repeat the numerical analysis for 
the 7r/5-right triangle billiard. In this case the length 
scale L max , for the saturation of the number variance, is 
about 2.6 times smaller than for the barrier billiard, i.e. 

L max = ^2^/(8 sin 2a) a = n/5. (16) 
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(Tl) (T2) (T3) 

jV = 4 * 10 5 10 6 4 * 10' 
Ltot = 4 * 10 4 10 5 10 5 



TABLE II: Data sets for the right triangle billiard (a = 7r/5). 
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FIG. 2: The scaled number variance for the data sets in 
Tab.[H] Only a subset of the data points is plotted. The solid 
line gives the fit g(x) with parameters as given in Eq. 1171 . 

We calculated three data sets in different energy regions, 
as listed in Tab. |HJ The stepsize was again L st = 10. 
With all data sets, we fit the scaled number variance 
g(x) using the parametrization (|14(l . excluding as before 
all data points with x > 0.2. The resulting estimates are: 

a = 0.5385(57) ai = 0.958(37) 

a 2 = 0.4154(54) a 3 = 46.2(5.4) (17) 

with a reduced chi-square value of xlt// ~ 0.20 
(/ = 61). The numerical results for g(x), and the fit 
with Eq. (|14(l . are plotted in Fig. [21 For the error bars, 
we have checked, that the same approximation holds 



as in the case of the barrier billiard. The absolute 
errors are smaller here, because the data sets are taken 
at lower energy ranges. As in Fig. Q we find two 
different linear regimes, with a transition which occurs 
somewhat earlier, at i « 0.04. The overall change in 
the graph is less dramatic because the slope at x < 0.04 
is much smaller in magnitude. Though the agreement 
of the extrapolated value for g(0) with the theoretical 
expectation ^ = 5/9 ^| is not perfect, our estimate is 
quite close to it. Note that the error estimates for the 
parameters clq, . . . , 03 are also based on statistical data. 
In particular in the case of do, there are only relatively 
few relevant data points. This leads to rather large 
uncertainties. In general the data points fluctuate much 
less than in the case of the barrier billiard, as can be 
read off from the smaller value for x| t //- 

We advanced the technique, recently proposed in 0, to 
solve the Schrodinger equation for perturbations of low 
rank. We realized that due to a special sum rule, one 
can obtain the level counting function with very little 
numerical effort. Systems involving short range interac- 
tions may be good candidates for future applications. In 
this letter, we considered two pseudointegrable billiards: 
the barrier billiard and the 7r/5-right triangle billiard. In 
spite of their apparent simplicity, the spectral statistics, 
a corner stone in quantum-classical correspondence, is 
only hardly understood. We performed extensive numer- 
ical calculations for the number variance E 2 (L) , L 1 
in energy regions up to M > 1.6 * 10 7 (for the barrier 
billiard) and Af > 4. * 10 6 (for the right triangle billiard). 
With the help of the scaled number variance, we obtained 
very precise estimates for the level compressibility, which 
largely confirm previously obtained analytical results. 
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